El siguiente documento corresponde al primer avance del trabajo de investigación del curso Visualización de datos y análisis estadístico del Posgrado Analista de Datos de la Universidad Cenfotec.
A continuación, se exponen el tema a tratar en esta investigación, los objetivos planteados, asi como una muestra del dataset que se empleara en dicha investigación.
“Estadisticas y Modelado de la distribución biogeográfica de Pumas y Jaguares de un sector del continente americano”
https://www.inaturalist.org/
https://biodiversityinformatics.amnh.org/
library(dplyr)
library(tibble)
library(tidyr)
library(sp)
library(DT)
library(leaflet)
library(leaflet.extras)
library(raster)
library(rasterVis)
library(rgdal)
library(ggplot2)
library(plotly)
library(LaCroixColoR)
library(sf)
library(here)
library(devtools)
library(proto)
library(highcharter)
library(rgeos)
library(velox)
library(usdm)
library(gtools)
library(corrplot)
library(Hmisc)
library(dismo)
library(ENMeval)
En esta sección se realizan una serie de procesos que conllevan al desarrollo de las bases de datos finales
Carga de la base de datos y de las variables que la conforman
df <- read.csv("E:/Analista de datos/DAT004/Proyecto Final/Datasets/samples/wildcats_df.csv", stringsAsFactors=F, sep=",", na.strings=c("","NA"), encoding="UTF-8")
colnames(df)
## [1] "X.U.FEFF.common_name" "scientific_name" "latitude"
## [4] "longitude" "created_at" "updated_at"
## [7] "place_country_name"
Cambio de nombre de las columnas
df <- df[,c(1:5,7)]
names(df)[1] <- "Nombre_Comun"
names(df)[2] <- "Nombre_Cientifico"
names(df)[3] <- "Latitud"
names(df)[4] <- "Longitud"
names(df)[5] <- "Fecha_Registro"
names(df)[6] <- "Pais"
colnames(df)
## [1] "Nombre_Comun" "Nombre_Cientifico" "Latitud"
## [4] "Longitud" "Fecha_Registro" "Pais"
Se imprime un resumen de datos por especie de modo que se pueda obtener un conteo inicial
df%>%
dplyr::group_by(Nombre_Comun)%>%
dplyr::summarise(length(Nombre_Comun))
Se realiza una correccion de datos en la columna de “Nombre_Comun” de forma que se estandaricen algunos registros que figuran como duplicados
ndf <- df%>%
mutate(Nombre_Comun=ifelse(Nombre_Comun=="Puma de América del Norte", "Puma",Nombre_Comun))%>%
mutate(Nombre_Comun=ifelse(Nombre_Comun=="Puma de América del Sur", "Puma",Nombre_Comun))
Nuevamente se imprime un resumen final de datos por especie
count_sp <- ndf%>%
group_by(Nombre_Comun)%>%
summarise(Count=length(Nombre_Comun))
plot_ly(data=count_sp, x = ~Nombre_Comun, y = ~Count, type = 'bar', text = ~Count, textposition = 'auto',
marker = list(color = c('rgba(204,204,204,1)', 'rgba(222,45,38,0.8)'),
line = list(color = c('rgba(204,204,204,1)', 'rgba(222,45,38,0.8)'),
width = 1.5))) %>%
layout(title = "Cantidad de especies registrados en el dataframe inicial",
xaxis = list(title = "Especies"),
yaxis = list(title = "Cantidad"))
cf <- colorFactor(c("#ffa500", "#13ED3F"), domain=ndf$Nombre_Comun)
m <- leaflet(df) %>% addTiles(group = "OSM (default)") %>%
addProviderTiles(providers$Stamen.Toner, group = "Toner") %>%
addProviderTiles(providers$Stamen.TonerLite, group = "Toner Lite") %>%
setView(-63.54105072, -8.88123188, zoom = 2) %>%
addCircles(~Longitud, ~Latitud, popup=ndf$Nombre_Comun, group = ndf$Nombre_Comun ,weight = 3, radius=4,
color=cf(ndf$Nombre_Comun), stroke = TRUE, fillOpacity = 0.5) %>%
addLegend("bottomright", colors= c("#ffa500", "#13ED3F"), labels=c("Puma", "Jaguar"), title="Especies") %>%
addLayersControl(
baseGroups = c("OSM (default)", "Toner", "Toner Lite"),
overlayGroups = ndf$Nombre_Comun,
options = layersControlOptions(collapsed = TRUE))%>%
suspendScroll()
m
El area de estudio con la que se pretende trabajar corresponde a Centro y Sudamrica, asi como el Caribe.
Algunas de las coberturas que se utilizan para conocer el habitat de las especies son variables climáticas que se derivan de los datos proporcionados por el Panel Intergubernamental sobre Cambio Climático y se produjeron utilizando Interpolación de basada en lecturas tomadas en estaciones meteorológicas de todo el mundo desde 1961 hasta 1990.
Los coberturas corresponden a:
*cld6190_ann : cobertura de nubes, anual
*dtr6190_ann : rango de temperatura diurna, anual
*frs6190_ann : frecuencia de heladas, anual
*pre6190_ann : precipitacion, anual
*pre6190_I1 : precipitacion, enero
*pre6190_I4 : precipitacion, abril
*pre6190_I7 : precipitacion, julio
*pre6190_I10 : precipitacion, octubre
*tmn6190_ann : temperatura promedio, anual
*tmp6190_ann : temperatura minima, anual
*tmx6190_ann : temperatura maxima, anual
*vap6190_ann : presion de vapor, anual
r_files <- list.files(here::here("Datasets", "coverages"),
full.names = T)
r_files
## [1] "E:/Analista de datos/DAT004/Proyecto Final/Datasets/coverages/cld6190_ann.asc"
## [2] "E:/Analista de datos/DAT004/Proyecto Final/Datasets/coverages/dtr6190_ann.asc"
## [3] "E:/Analista de datos/DAT004/Proyecto Final/Datasets/coverages/ecoreg.asc"
## [4] "E:/Analista de datos/DAT004/Proyecto Final/Datasets/coverages/frs6190_ann.asc"
## [5] "E:/Analista de datos/DAT004/Proyecto Final/Datasets/coverages/h_dem.asc"
## [6] "E:/Analista de datos/DAT004/Proyecto Final/Datasets/coverages/pre6190_ann.asc"
## [7] "E:/Analista de datos/DAT004/Proyecto Final/Datasets/coverages/pre6190_l1.asc"
## [8] "E:/Analista de datos/DAT004/Proyecto Final/Datasets/coverages/pre6190_l10.asc"
## [9] "E:/Analista de datos/DAT004/Proyecto Final/Datasets/coverages/pre6190_l4.asc"
## [10] "E:/Analista de datos/DAT004/Proyecto Final/Datasets/coverages/pre6190_l7.asc"
## [11] "E:/Analista de datos/DAT004/Proyecto Final/Datasets/coverages/tmn6190_ann.asc"
## [12] "E:/Analista de datos/DAT004/Proyecto Final/Datasets/coverages/tmp6190_ann.asc"
## [13] "E:/Analista de datos/DAT004/Proyecto Final/Datasets/coverages/tmx6190_ann.asc"
## [14] "E:/Analista de datos/DAT004/Proyecto Final/Datasets/coverages/vap6190_ann.asc"
st <- raster::stack(r_files)
raster::plot(st)
ecoreg = raster(r_files[3])
raster::plot(ecoreg, main="Ecoregiones del area de estudio")
raster::hist(ecoreg, main="Histograma de las Ecoregiones del area de estudio", col="green")
raster::pairs(st[[c(1,6,11,14)]])
A continuacion, se carga el un archivo espacial de tipo vectorial correspondiente a los paises que conforman el area de estudio.
area_estudio <- sf::st_read(dsn="E:/Analista de datos/DAT004/Proyecto Final/Datasets/samples/map.geojson", layer="map", quiet=TRUE)
map.area <- plot_geo(area_estudio, name=~pais, color = I("#89C5DA"))%>% hide_legend()
map.area
De previo se realiza una selección de la base de datos de las especies respecto del área de estudio
coords <- data.frame(x=ndf$Longitud, y=ndf$Latitud)
prj <- CRS("+proj=longlat +datum=WGS84 +no_defs")
points <- SpatialPointsDataFrame(coords, ndf, proj4string = prj)
wildcats <- crop(points,ecoreg)
wildcats <- data.frame(wildcats)
wildcats <- wildcats[,1:6]
Para luego visualizar los puntos de observacion de las especies a estudiar
cf2 <- colorFactor(c("#ffa500", "#13ED3F"), domain=wildcats$Nombre_Comun)
m2 <- leaflet(wildcats) %>% addTiles(group = "OSM (default)") %>%
addProviderTiles(providers$Stamen.Toner, group = "Toner") %>%
addProviderTiles(providers$Stamen.TonerLite, group = "Toner Lite") %>%
setView(-63.54105072, -8.88123188, zoom = 2) %>%
addCircles(~Longitud, ~Latitud, popup=wildcats$Nombre_Comun, group = wildcats$Nombre_Comun ,weight = 3, radius=4,
color=cf2(wildcats$Nombre_Comun), stroke = TRUE, fillOpacity = 0.5) %>%
addLegend("bottomright", colors= c("#ffa500", "#13ED3F"), labels=c("Puma", "Jaguar"), title="Especies") %>%
addLayersControl(
baseGroups = c("OSM (default)", "Toner", "Toner Lite"),
overlayGroups = wildcats$Nombre_Comun,
options = layersControlOptions(collapsed = TRUE))%>%
suspendScroll()
m2
Primero se realiza una selección final de la base de datos e inclusion de nuevas variables
wildcats$Hora_Registro <- strftime(wildcats$Fecha_Registro, format="%H") #col7
wildcats$Fecha_Registro <- as.POSIXct(wildcats$Fecha_Registro, tz="", format="%Y-%m-%d") #col5
wildcats$Y.M.R <- strftime(wildcats$Fecha_Registro, format="%Y-%m") #col8
wildcats$M.R <- format(wildcats$Fecha_Registro, "%m") #col9
wildcats$Y.R <- format(wildcats$Fecha_Registro, "%Y") #col10
wildcats <- wildcats[, colnames(wildcats)[c(1:4,6,5,8,10,9,7)]] #Reorganizacion de las columnas
Se procede a visualizar la base de datos
datatable(wildcats, extensions = c('Buttons','FixedColumns', 'RowReorder'),
filter = list(position = 'top', clear = FALSE),
options = list(dom = 'Bfrtip',
buttons = c('copy', 'csv', 'excel', 'pdf', 'print'), dom = 't',
scrollX = TRUE,
fixedColumns = TRUE,
rowReorder = TRUE, order = list(c(0 , 'asc'))))
En este primer apartado se realizan analisis sobre estadísticas espacio-temporales respecto de los pumas y jaguares de un sector del continente americano
plot_ly(x= wildcats$Pais,split=wildcats$Nombre_Comun, color= I(cf2(wildcats$Nombre_Comun)),
frame=wildcats$M.R, type = "histogram", alpha = 0.6)%>%
layout(barmode = "overlay")
## Warning: Ignoring 3 observations
c.jaguars <- wildcats%>%
dplyr::filter(Nombre_Comun == "Jaguar")%>%
dplyr::group_by(Pais)%>%
dplyr::summarise(Cant=n())
c.pumas <- wildcats%>%
dplyr::filter(Nombre_Comun == "Puma")%>%
dplyr::group_by(Pais)%>%
dplyr::summarise(Cant=n())
data(worldgeojson, package = "highcharter")
highchart() %>%
hc_title(text = "Cantidad de Jaguares") %>%
hc_subtitle(text = "por pais") %>%
hc_add_series_map(map = worldgeojson, df = c.jaguars, name = "cantidad de individuos",
value = "Cant", joinBy = c("name", "Pais"), allAreas=FALSE, animation=TRUE)%>%
hc_colorAxis(min= 1,
type='logarithmic',
stops = color_stops())%>%
hc_tooltip(useHTML = TRUE, headerFormat = "",
pointFormat = "Pais: {point.name}/ Cantidad de especies: {point.Cant}")%>%
hc_legend(layout= 'horizontal',
verticalAlign= 'top')%>%
hc_mapNavigation(enabled = TRUE)
data(worldgeojson, package = "highcharter")
highchart() %>%
hc_title(text = "Cantidad de Pumas") %>%
hc_subtitle(text = "por pais") %>%
hc_add_series_map(map = worldgeojson, df = c.pumas, name = "cantidad de individuos",
value = "Cant", joinBy = c("name", "Pais"), allAreas= FALSE, animation=TRUE)%>%
hc_colorAxis(min= 1,
type='logarithmic',
stops = color_stops())%>%
hc_tooltip(useHTML = TRUE, headerFormat = "",
pointFormat = "Pais: {point.name}/ Cantidad de especies: {point.Cant}")%>%
hc_legend(layout= 'horizontal',
verticalAlign= 'top')%>%
hc_mapNavigation(enabled = TRUE)
c.jaguars <- wildcats%>%
dplyr::filter(Nombre_Comun == "Jaguar")%>%
dplyr::group_by(Fecha_Registro)%>%
dplyr::summarise(Cant=n())
c.pumas <- wildcats%>%
dplyr::filter(Nombre_Comun == "Puma")%>%
dplyr::group_by(Fecha_Registro)%>%
dplyr::summarise(Cant=n())
ct.jaguars <- ts(data=c.jaguars, start = c(2011,10), end = c(2019,12), frequency = 12)
ct.pumas <- ts(data=c.pumas, start = c(2011,10), end = c(2019,12), frequency = 12)
highchart(type = "stock") %>%
hc_title(text = "Cantidad de especies") %>%
hc_subtitle(text = "en la linea de tiempo de observacion") %>%
hc_add_series(ct.jaguars[,2], type = "line", pointInterval= 30*24*3600*1000, name="Cant de jaguares observados a la fecha", color = "#ffa500")%>%
hc_add_series(ct.pumas[,2], type = "line", pointInterval= 30*24*3600*1000, name="Cant de pumas observados a a fecha", color = "#13ED3F")%>%
hc_xAxis(title= list(text='Fechas'), type='datetime')%>%
hc_yAxis(title= list(text='Cantidad de especies observadas'), opposite = FALSE)
decomp_jaguars <- decompose(ct.jaguars[,2])
decomp_pumas <- decompose(ct.pumas[,2])
plot(decomp_jaguars)
par(adj = 1)
title(sub = "Jaguares", cex.sub = 1.5, font.sub = 3, col.sub = "#ffa500")
plot(decomp_pumas, sub="Pumas")
par(adj = 1)
title(sub = "Pumas", cex.sub = 1.5, font.sub = 3, col.sub = "#13ED3F")
e.m <- ggplot(wildcats, aes(x=Longitud, y=Latitud, frame=M.R, text = paste("Pais:", Pais))) +
borders("world", xlim = c(-120, -30), ylim = c(-50, 25))+
geom_point(data=filter(wildcats), shape = 21, colour = cf2(wildcats$Nombre_Comun), size = 1)+
facet_grid(cols = vars(Nombre_Comun)) +
labs(title = "Distribucion de las especies observadas segun mes")
ggplotly(e.m)%>%
animation_opts(
2000, easing = "elastic", redraw = FALSE
)%>%
animation_slider(
currentvalue = list(prefix = "Mes ", font = list(color="red"))
)
pd <- qplot(data = wildcats, x=Longitud, y=Latitud) +
stat_ellipse(geom = "polygon", alpha = 1/2, aes(fill = Nombre_Comun, frame=M.R))+
borders("world", xlim = c(-120, -30), ylim = c(-50, 25))+
geom_point(data=filter(wildcats), aes(frame=M.R), shape = 21, fill = cf2(wildcats$Nombre_Comun), size = 2)+
facet_grid(cols = vars(Nombre_Comun)) +
labs(title = "Patron de Distribucion de las Especies")
ggplotly(pd)%>%
animation_opts(
2000, easing = "elastic", redraw = FALSE
)%>%
animation_slider(
currentvalue = list(prefix = "Mes ", font = list(color="red"))
)
dn <- ggplot(wildcats, aes(x=Longitud, y=Latitud)) +
borders("world", regions=".", xlim = c(-120, -30), ylim = c(-50, 25))+
stat_density2d(aes(fill = stat(level), frame=wildcats$M.R), geom="polygon", alpha=0.3) +
geom_point(data=filter(wildcats), aes(frame=M.R), shape = 21, fill = cf2(wildcats$Nombre_Comun), size = 2, alpha=0.3)+
facet_grid(cols = vars(Nombre_Comun)) +
scale_fill_viridis_c(option = "viridis") +
theme(legend.position = "magma") +
labs(title = "Densidad segun distribucion por especie")
ggplotly(dn)%>%
animation_opts(
4000, easing = "elastic", redraw = FALSE
)%>%
animation_slider(
currentvalue = list(prefix = "Mes ", font = list(color="red"))
)
ecoreg_spec <- raster::extract(ecoreg, wildcats[,4:3]) %>%
cbind(wildcats, .) %>%
as.data.frame()
ecoreg_spec <- na.omit(ecoreg_spec)
ecoreg_spec$Ecoreg <- as.character(ecoreg_spec$.)
ggplot(ecoreg_spec, aes(x=Longitud, y=Latitud, color = ecoreg_spec$Ecoreg)) +
borders("world", xlim = c(-120, -30), ylim = c(-50, 25))+
geom_point(data=filter(ecoreg_spec), shape = 21, size = 1)+
scale_colour_manual(values = lacroix_palette("PassionFruit", n=14))+
facet_grid(cols = vars(Nombre_Comun)) +
labs(title = "Especies observadas segun ecoregion en la que se encuentran")
ggplot(ecoreg_spec, aes(y=Ecoreg, x=Nombre_Comun, colour = Ecoreg)) +
geom_count(alpha=0.5) +
scale_colour_manual(values = lacroix_palette("PassionFruit", n=14))+
labs(title = "Cantidad de especies observadas por ecoregion",
x = "Especie",
y = "Ecoregiones",
size = "")
Proceso de union de los puntos de observacion de las especies y los raster de las variables ambientales
#Se extraen temas relevantes de las variables ambientales
select_r = dropLayer(st,3)
jp_swd <- raster::extract(select_r, wildcats[,4:3]) %>%
cbind(wildcats, .) %>%
as.data.frame()
#Se elimina las filas con valores nulos
jp_swd <- na.omit(jp_swd)
Ploteo de la correlacion de variables
mt <- cor(jp_swd[,11:23])
corrplot(mt, method = 'circle', type = 'upper')
Esto se hace a traves del analisis en regresion lineal VIF, el cual es el factor de inflación de la varianza. Este cuantifica la intensidad de la multicolinealidad en un análisis de regresión normal de mínimos cuadrados. Proporciona un índice que mide hasta qué punto la varianza (el cuadrado de la desviación estándar estimada) de un coeficiente de regresión estimado se incrementa a causa de la colinealidad.
# Analisis VIF
vif.res <- vif(x = jp_swd[,11:ncol(jp_swd)])
vif.step <- vifstep(x = jp_swd[,11:ncol(jp_swd)], th = 10)#Los factores VIF mayores que 10 implican problemas graves de multicolinealidad
vrs <- vif.step@results$Variables %>% as.character()
vif.step
## 4 variables from the 13 input variables have collinearity problem:
##
## pre6190_ann tmp6190_ann tmn6190_ann tmx6190_ann
##
## After excluding the collinear variables, the linear correlation coefficients ranges between:
## min correlation ( pre6190_l10 ~ pre6190_l1 ): -0.009063609
## max correlation ( vap6190_ann ~ frs6190_ann ): -0.8508242
##
## ---------- VIFs of the remained variables --------
## Variables VIF
## 1 cld6190_ann 3.025670
## 2 dtr6190_ann 1.678499
## 3 frs6190_ann 4.102199
## 4 h_dem 2.447628
## 5 pre6190_l1 1.755003
## 6 pre6190_l10 5.388244
## 7 pre6190_l4 4.074227
## 8 pre6190_l7 2.896292
## 9 vap6190_ann 5.176151
vrs
## [1] "cld6190_ann" "dtr6190_ann" "frs6190_ann" "h_dem" "pre6190_l1"
## [6] "pre6190_l10" "pre6190_l4" "pre6190_l7" "vap6190_ann"
jp_fnl <- jp_swd %>%
as_tibble %>%
dplyr::select(Nombre_Cientifico:Longitud, vrs)
head(jp_fnl)